In the example above, we let interc_coeff='k' be the coefficient for all terms. This will be replaced with a ConstantKernel whose parameter value will be optimized. Also, interc_coeff here represents the independednt additional error that can be added to the kernel, typically WhiteNoise() + RBF().
The search for the optimal kernel polynomial is typically a non-convex optimization problem with many local minima/optimal kernels. In fact, a large number of models which can perfectly fit a finite dataset can be found. This can be limited by reducing the kernel space to search, based on existing knowledge of the process that generated the data.
We will use the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) algorithm, which is an evolutionary strategy (stochastic, derivative-free) that is available in Optuna as optuna.samplers.CmaEsSampler. Note that this sampler does not support categorical distributions. For these situations, we would opt for the TPESampler (Tree-structured Parzen Estimator). The latter optimization approach is described in details in Bergstra, James S., et al. “Algorithms for hyper-parameter optimization.” Advances in Neural Information Processing Systems. 2011 and Bergstra, James, Daniel Yamins, and David Cox. “Making a science of model search: Hyperparameter optimization in hundreds of dimensions for vision architectures.” Proceedings of The 30th International Conference on Machine Learning. 2013.
The search is set to explore 1000 trials with an early stopping if no improvement is found in 200 trials. The objective function builds and trains a modeling pipeline that preprocesses the data (scaling) and fits a Gaussian process regression model using the kernel of the trial. For smaller datasets (50 or fewer data rows), minimizing the Bayesian Information Criterion (BIC) drives the search for the optimal polynomial kernel. BIC is defined as
$BIC = k\log(n) - 2\log(\hat{L})$
where $\hat{L}$ is the estimated log-likelihood of the Gaussian process regression model, $n$ is the number of data rows, $k$ is the number of parameters of the kernel function. The hyperparameter values of the kernel are optimized when fitting the Gaussian process regressor by maximizing the log-marginal-likelihood (LML) as described here.
For larger datasets, adopting the approach described above would lead to impracticable computational times. Therefore, we stratify the dataset based on the target variable using KBinsDiscretizer. The target values are binned so that all bins have the same number of data (quantile strategy). The stratification happens by utilizing the bin IDs. At each trial, 50 data rows are sampled using this stratification of the target and are used to train the Gaussian process regressor while optimizing its kernel hyperparameter values. Different data samples are utilized at each trial. In the case of a kernel already tried previously in the search, the kernel is initialized with the hyperparameter values found earlier. The whole dataset is then used to only fit the model without further optimizing its kernel hyperparameter values. The search of the optimal polynomial kernel is driven by minimizing BIC calculated using the LML from this fit. A fundamental challenge with Gaussian processes is scalability. The covariance matrix inversion complexity increases quadratically with the size of the data.
To test this approach, we let Optuna find an optimal polynomial with no more than 2 base kernels in product terms and no more than 3 total terms plus the additional noise term. This will allows us to see how the result compares to the kernel used in the study described in the previous section.
Running multiple times may result into slightly different kernel polynomials since this optimization problem is not convex. In fact, we should always look at the top trials to make a decision on the kernel to adopt based on our intuition of the problem.
